Time-Dependent Response Calculations of Nuclear Resonances 
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A new alternate method for evaluating linear response theory is formally developed, and results 
are presented. This method involves the time-evolution of the system using TDHF and is constructed 
directly on top of a static Hartree-Fock calculation. By Fourier transforming the time-dependent 
result the response function and the total probability amplitude are extracted. This method allows 
for a coherent description of static properties of nuclei, such as binding energies and deformations, 
while also providing a method for calculating collective modes and reaction rates. A full 3-D 
Cartesian Basis-Spline collocation representation is used with several Skyrme interactions. Sample 
results are presented for the giant multipole resonances of ^®0, *'^Ca, and ^^S and compared to other 
calculations. 

PACS numbers: 21.60.-n,21.60. Jz 



I. INTRODUCTION 

Typically, linear response theory equations are derived 
by adding a specific time-dependent perturbing function 
to the Hamiltonian, which is usually harmonic in time, 
resulting in a set of RPA-like equations. These equations 
are then solved, for a given energy, using several methods 
(see for example Refs. tL2e^M) to give the response of 
the system to a specific collective excitation mode. 

In this paper, a new alternate method is presented to 
calculate response theory. A specific time-dependent per- 
turbing external piece is added to the static Hamiltonian 
to give a time-dependent total Hamiltonian, Htot{t)- A 
static Hartree-Fock solution is then time-evolved using 
this Htot{t) in a time-dependent Hartree-Fock (TDHF) 
calculation. The time-dependent result is then Fourier 
transformed to give the response of the system for all 
energies. In this scheme one recovers both the response 
spectrum as well as the total transition probability am- 
plitude corresponding to a given specific collective mode. 
Similar analyses of the long-time evolution of TDHF 
equations to study collective vibrations have been uti- 
lized in the past 0, H 0, S 0, as well as extensions 
to study the damping of giant resonances [lillll. The 
main advantage of this approach is that the dynami- 
cal response calculation is constructed directly on top 
of a static Hartree-Fock calculation and hence the static 
and dynamical calculations are calculated using the same 
Hamiltonian description. Therefore there is a complete 
consistency between the static ground state of the system 
and the response calculations. One can then provide a 
coherent description of static properties of nuclei and of 
dynamical properties. This is important for example in 
/3— decay calculations of exotic nuclei, where reliable pre- 
dictions are verysensitive to the deformation properties 
of the nucleus Hence in this formalism consistent 

predictions of both the deformation and reaction rate 
properties are possible. 

The static and dynamical Hartree-Fock calculations 
are performed usin g a 3-D Cartesian Basis-Spline col- 
location expansion |l3l [T3 | . Properties of Basis-Splines 



include the practicality associated with coordinate-space 
lattice grids, while also providing accurate representa- 
tions of the gradient operator and a good description of 
the continuum. 

In Section II, the time-dependent evaluation of the re- 
sponse theory is discussed and shown to give the total 
transition probability. Numerical details and sample re- 
sults are presented in Section III. 



II. TIME-DEPENDENT RESPONSE THEORY 

The response equations can be derived from a specific 
time-dependent perturbation functional of the TDHF 
equations . To begin the proof a solution to the static 
Schrodinger equation is written as follows: 



i/|Vs(0)) = i?|Vs(0)). 



(1) 



A time-dependent perturbing function is added to the 
static Hamiltonian: 



Htot — H + Hex {t) ■ 
The external piece is defined as: 

(fx n{x, t)F{x) 



fit), 



(2) 



(3) 



where n{x,t) is the number density operator and F(x) 
corresponds to a one-body operator to excite a particular 
collective mode. The functions, F{x) and f{t) will be 
chosen later. 

At some time t — to the external piece of the Hamil- 
tonian is turned on where \ipsit)) is the solution to the 
time-dependent Schrodinger equation: 



(4) 



Here the subscript "s" refers to the Schrodinger picture. 
A solution of the following form is then constructed: 



(5) 
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where for t < to, A{t) = 1. Using Eq. the function 
A{t) can be shown to be a solution to: 

zh^Ait) = HUt)A{t), (6) 



where the superscript "7" refers to the interaction pic- 
ture, which reduces to the Heisenberg picture when 
Hex{t) = 0. The solution to Eq. © can be written itcr- 
atively as: 



A(t) = l-- / dt'Hi,{t') + ..., (7) 



where the state vector is then given by: 



\^,{t)) = e-^"'/'-^\MO)) - ^e-^"'/" j'dt'HUt')\MO)) + - ■ (8) 

The expectation value of any operator, 0{t), is equal to 

{Mt)\Os{t)\Mt)) = (^-.(0)107(01^,(0))+ (9) 

+(^^,(0)1 J ^ dt' [HUt'),di{t)] mo)) + ... . 

The linear approximation is made such that terms beyond first order in are neglected. The first term in the 
expansion is trivially the unperturbed expectation value of the operator in the Schrodinger picture. If we choose the 
operator 0{t) to be the number density operator, then using Eq. (PJ, the fiuctuation in the density can be defined as: 

S{n{x,t)) = {Mt)\ns{t)\Mt)) - {MO)\nsmMO)) 

= iMO)\^J^ dt' J d''x'F{x')f{t')[ni{x',t'),ni{x,t)]\MO)) ■ (10) 
The retarded density correlation function is defined as: 

zD^ix, t; x', t') = e{t - ^,^ {M[nH{xlnH{x'mo) ^ ^^^^ 

{mm) 

where Uh = tlh — {nn) is the deviation of the number operator in the Heisenberg picture. The density fiuctuation 
can be written as: 

S{n{x,t))^^J dt' J d^x'D"'{x,t;x',t')F{x')f{t'). (12) 
Using the Fourier representation of 9{t — t'), the Fourier transform of the density correlation function is: 



/oc 
d{t-t')e"^^'-''hD^{x,t; x',t') 
-OO 



^ \p J {'4'o\nsix)\'iJjn) {ijM\ns{x')\'ilJo) {ipo\nsix')\'ipn){'ipn\nsix)\'ipo) I 
~ n\ io-^^+zr^ ^ + J' 

where {ipn) represents the full spectrum of the excited many-body states of H . The Fourier transform of the density 
fiuctuation then becomes: 

/OO 
dte"^'5{n{x,t)) (14) 
-OO 



i J d^x'D"{x,x';Lu)Fix')fico), 



where 



fico) = / dt'e^-'' f{t') . (15) 
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The linear response structure function, S{lu) is derived to be: 

f{uj)S{uj) = / d^xS{F''{x)n{x,Lu)) 



\ f d^x f d^x'F\x)D^{x,x';uj)F{x')f{uj). 



(16) 

Combining Eqs. H13() and H16|) the imaginary part of the structure function then gives the total transition probability 
associated with F{x) 



Im [S{uj)] = -lY. J d''x'{i>n\nAx')\i'o)F{x': 



6 uj 



En — Eq 



En > Eo- 



(17) 



Note that this quantity is negative definite; this feature can be used as a measure of the convergence of the solution. 



At this point instead of using the standard route of let- 
ting f{t) — > to recover the linear response equations, we 
choose an alternate technique to calculate the response. 
In this case we evolve the system in time and then Fourier 
transform the result, where Hex{t) is a perturbing func- 
tion. We choose f{t) to be a Gaussian of the following 
form, 




t > to 



(18) 



where e is some small number (~ 10~^), chosen such that 
we are in the linear regime. The parameter, a, is set to 
be 1.0 c^/fm^, which allows for a reasonable perturbation 
of collective energies up to « 150 MeV. 

In practice, our numerical calculations proceed as fol- 
lows: First, we generate highly accurate static HF wave 
functions on the 3-D lattice. Then the external time- 
dependent perturbation, Eq. (PJ, is set up by choosing 
a particular form for F{x), and using Eq. ((TH|) for f{t). 
Next, we solve the TDHF equations utilizing the time- 
evolution operator 



U{t,to)=T 



-j- r dt'Htotxt') 



(19) 



where T[ ... ] denotes time-ordering. Using infinitesimal 
time increments, the time-evolution operator is approxi- 
mated by 



U{tn+i,t„) = e « 



J'"+^ dt'Htotit') 



(20) 



TV 



fc=i 



kl 



where the quantity H^^^^ is evaluated by repeated oper- 
ations of Hfot upon the wave functions. From the nu- 
merical solution of the TDHF equations, we obtain the 
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FIG. 1: The fluctuations in the isoscalar axial quadrupole 
moment as a function of time are shown for ^^O, ^^S, and 
^''Ca using the SkM* interaction. The size of the time step 
is 0.4 fm/c. A grid of 24'^ with a Cartesian box dimensioned 
(-12,-hl2 fm)=' is used. 



density fluctuation as a function of time, Eq. H12|l . Af- 
ter Fourier transforming this quantity and folding with 
the function F{x), we obtain f{uj)S{uj) in Eq. from 
which the linear response structure function of the sys- 
tem is extracted. 



III. NUMERICAL DETAILS AND RESULTS 

The static and time-dependent Hartree-Fock calcula- 
tions are performed using a collocation spline basis in 
a three-dimensional lattice configuration. Basis-Splines 
allow for the use of a lattice grid representation of the 
nucleus, which is much easier to use than alternate ba- 
sis techniques, such as multi-dimensional harmonic os- 
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dilators. Also, for studies of exotic nuclei, because of 
weak binding, the density distributions tend to extend 
to large distances and hence one finds a large sensitiv- 
ity to a harmonic oscillator basis due to the unphysical 
description of continuum states, while for a lattice grid 
representation one needs to simply increase the size of 
the box. Traditional grid representations typically use 
finite difference techniques to represent the gradient op- 
erator. Collocation Basis-Splines allow for the gradient 
operator to be represented by its action upon a basis 
function in a matrix form. Thus the collocation method 
gives a much more accurate representation of the gradi- 
ent, while maintaining the convenience of a lattice grid 
and hence provides a much more accurate calculation in 
the end 
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FIG. 2: The imaginary part of the response function corre- 
sponding to the isoscalar quadrupole moment is shown for ^^O 
using three different Skyrme force parametrizations, where 
the SkM*, Sgll and Skll versions are shown in the upper, 
middle and lower panels, respectively. The experimentally 
measured giant quadrupole resonance is around 20.7 MeV. 

We have performed the usual tests for assuring con- 
vergence with respect to the numerical box size and the 
number of collocation points. Typically a 7th order B- 
spline is used in a (12 fm)'^ box with 20"^ — 24^^ grid points. 
The calculation can be performed with or without as- 
suming time-reversal symmetry. The collective linear re- 
sponse may involve particle-hole interactions which are 
spin-dependent and not time-reversal symmetric. There- 
fore, for the correct collective content to be included one 
should not impose time-reversal symmetry in the linear 
response calculations In the results presented in this 
paper time-reversal symmetry is not imposed. In a com- 
parison it is found that imposing time-reversal symmetry 
causes small shifts in the position of the collective modes 
on the order of ~ 0.3 MeV. 

Calculations of isovector dipole; isovector and isoscalar 
octupole; and isoscalar quadrupole collective modes are 
performed for ^^O, '^^S, and ^"Ca using several parame- 



terizations of the Skyrme interaction. Here the param- 
eterizations known as Skll [3, SkM* |l3 and Sgll [i| 
are used for comparisons. The exponent of the density 
in the density-dependent term in Skll is a = 1, while 
for SkM* and Sgll this exponent is a = i. This causes 
the Skll force to produce a rather large nuclear matter 
incompressibility, while the SkM* and Sgll forces pro- 
duce more realistic compression properties. Also the 
SkM* and Sgll forces allow for more stable static Hartree- 
Fock and TDHF computations. The static Hartree-Fock 
calculations converge more easily and rapidly for SkM* 
and Sgll than for Ski, Skll or Skill. A time step of 
At — 0.4 fm/c is used for the calculation. It is found 
that one can perform the time evolution for up to 32768 
time steps without appreciable dissipation. The results 
shown here use 16384 time steps for a maximum time of 
w 6550 fm/c. 
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FIG. 3: FFT results for ^''O using three different numbers of 
time steps. The case when the number of time steps equals 
an exact power of two (32768) is clearly the better converged 
result, showing almost no positive values. 

Reasonable results are obtained if the parameter e in 
Eq. ((TH|l is chosen to fall in the range 2.0 x 10"^ < e < 
2 X 10"^. By varying the value of e, the amplitude of 
the time-dependent density fluctuation then scales pro- 
portionally to e, thus indicating that we are well within 
the linear regime of the theory. 

The linear response calculations require well converged 
initial static HF solutions. To test for the convergence of 
the static HF calculation the energy fluctuation, which is 
the variance of iJ, is minimized. The energy fluctuation 
is defined as 



(21) 



which measures how close the wave functions are to being 
eigenstates of H. This measure of convergence is very 
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sensitive to the eigensolutions and is independent of the 
iteration step size. For ^^O it was found that static HF 
solutions with energy fluctuation less than about 1.0 x 
10~^ provided adequate starting points for the dynamic 
calculation, although the smaller the energy fluctuation 
the better. 

The dynamic calculations involve using Eq. (|20(l to 
evolve the system. Since C/(t, t') is an unitary operator, 
the orthonormality of the system is preserved, therefore 
it is not necessary to re-orthogonalize the solutions af- 
ter every time-step. The stability of the calculation is 
checked by testing the preservation of the norm of each 
wave function. The number of terms in the expansion of 
the exponent in Eq. (|20|l is determined by requiring the 
norm to be preserved to a certain accuracy (typically to 
< 1.0 X 10-8 - 5.0 X 10-1°). 
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FIG. 4: The imaginary part of the response function cor- 
responding to the isoscalar quadrupole moment is shown for 
*°Ca using two different Skyrme force parametrizations, SkM* 
and Sgll. 



The time-dependent perturbing part of the Hamilto- 
nian is evaluated when the exponential term in Eq. I|18|) 
is greater than some small number, Scut- Since it is not 
difflcult to evaluate the action of the external part of the 
Hamiltonian on the wave function, Scut is chosen to be 
very small, (1.0 x 10^^'^). To allow the Fourier transform 
of f{t) to be evaluated easily, it is necessary to integrate t 
from — CX3 to c» and hence we would like the entire Gaus- 
sian of the perturbing function, f(t), to be included into 
the time evolution to the desired accuracy. The parame- 
ter to is therefore chosen such that the complete nonzero 
contribution of the time-dependent perturbation is in- 
cluded 
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(22) 



A. Quadrupole Excitation Modes 

For the study of the isoscalar quadrupole moment, 
the perturbing function F{x), introduced in Eq. JSJ, 
is chosen to be the mass quadrupole moment, Q20 = 
2z^ — {x^ + y^). It turns out that other even multipole 
modes are also excited at the same time (i.e. Q40, Qeo, 
...). One can therefore study the effect of the coupling 
between the different excitation modes. The same holds 
true for the odd multipoles. This is due to the nonlinear 
response effects present in the TDHF time evolution ^] . 

The quadrupole collective resonances are calculated for 
using three different Skyrme force parameterizations 
for comparisons. Two different size grids (22^^, 24'^) were 
used inside a (—12,-1-12 fm)'^ Cartesian box along with 
periodic boundary conditions. 
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FIG. 5: The imaginary part of the response function corre- 
sponding to the isoscalar quadrupole moment is shown for '^^S 
using the SkM* force. 

In Fig. ^the time-dependent evolution of the multipole 
moment defined as: 



(Q2o(t)) = / d^x6{n{x,t))Q2oix) , 



(23) 



is shown for ^^O, ^^S, and ^"Ca using the SkM* interac- 
tion. This figure illustrates the periodic character of the 
calculations with almost no damping. In this case the 
smallest oscillation is about 65 fm/c. 

A fast Fourier transform (EFT) is used to calculate 
the Fourier transform of {Q2o{t)) to give (Q2o(i^)) = 
/(cj)S'2o(t^). The time-dependent perturbation function, 
f{uj) can then be easily factored out using Eq. (|16ll . One 
can use the analytic expression for /(cj) or use a numer- 
ical FFT calculation, where the difference between the 
two methods ends up being negligible. 

In Fig. [3 the quadrupole responses for ^^O are shown 
for the three different Skyrme cases. The imaginary part 
of Eq. H23|) is Fourier transformed and divided by /(w). 
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Recalling Eq. H17|l . this quantity is derived to be a neg- 
ative definite quantity. The SkM* results for ^^O reflect 
this property very well in Fig. El We observe a sharp 
peak at about 20 MeV and a response which is almost 
purely negative. The side peaks, which are much less 
prominent, may not represent physical effects, but may 
be due to the construction of the continuum. In per- 
forming FFT transformations we find it to be important 
to have the number of points to be a power of two. In 
Fig. Elwe demonstrate this by plotting the response func- 
tion for three different choices for the total number of 
time steps which differ from each other by only 20 points 
each. As one can see the case for 32768 points results in 
almost no positive values, indicating good convergence. 
The three cases give similar pole structures, which are 
near the experimental isoscalar quadrupole giant reso- 
nance. The experimental peak is centered at an energy 
of 20.7 MeV with a width of about 7.5 ± 1 MeV 18]. The 
resonance calculated with the Sgll Skyrme parametriza- 
tion is closest to the experimental result, although SkM* 
is also fairly close. The Skll parametrization poorly rep- 
resents the data, but this force is known to give a bad rep- 
resentation of the collective nuclear properties. However, 
for the Skll parametrization we were able to compare our 
results for ^^O to the continuum RPA calculations of 
We find that the difference between the two calculations 
is less than 2% for the L=2, 3, and 4 modes. The differ- 
ence may be due to the omission of the spin-orbit term 
in the continuum RPA calculations. We find the closest 
agreement for the hexadecupole mode, where both cal- 
culations give a peak energy of 28.2 MeV. The widths 
of resonances are not accurately reproduced, because the 
continuum is included into the calculation in an approx- 
imate fashion and higher order correlations are missing 
in the TDHF approach. Since the calculations are per- 
formed in a box, the continuum is represented in terms of 
discrete pseudo-continuum states, whose density is sensi- 
tive to the size of the box. A larger box size is expected 
to better represent the continuum with a higher density 
of pseudo-continuum states. 

The property represented in Eq. (|17|l of being purely 
negative is not strictly reflected by the Sgll nor the Skll 
results. One can see that this property is approximately 
reflected, but it is clear the convergence of the solution 
for these two cases is not nearly as good as for the SkM* 
case. 

An alternate check of the calculation is the energy- 
weighted sum rule (EWSR). This provides a stringent 
test of the normalization and is derivable from the static 
Hamihonian. The SkM* resuh gives 92% of the EWSR, 
indicating excellent convergence. For Sgll and Skll the 
linear response results give 66% and 36% of the EWSR, 
respectively, indicating a lack of convergence for these 
results. 



The result for the Sgll calculation is particular puz- 
zling, since it is expected that the SkM* and Sgll inter- 
actions should be similar in behavior. The violations of 
the two convergence tests may be due to coupling to other 
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FIG. 6: The imaginary part of the response function corre- 
sponding to the isoscalar and isovector octupole, and isovector 
dipole are shown for ^"O using the SkM* force. 



collective modes, or may be due to numerical inconsisten- 
cies. We have investigated the Sgll result by varying the 
size of the box, while keeping the grid spacing the same. 
We have found that a shghtly larger box, (-14,-1-14) fm, 
seems to be closer to being purely negative. 

In Fig. 21 we show the quadrupole strength function 
for '^''Ca using two different Skyrme parameterizations. 
In this case SkM* and Sgll show a peak at about 17 
MeV. We also note the strength of the peak. In this case 
the resonance consumes all of the EWSR. In Fig. Owe 
show the same quantity for ■^^S using the Skyrme M* 
force. In this case we observe multiple structures in the 
giant quadrupole resonance; the energy splitting arises 
from the prolate quadrupole deformation of the nuclear 
ground state in this system. 

The resonant structure of the hexadecupole moment, 
(540, can also be calculated at the same time as the 
quadrupole moment. This corresponds to the coupling 
between the quadrupole and the hexadecupole collective 
modes. In this case if one were to include both Q2Q and 
Q40 into the external time-perturbing piece of the total 
Hamiltonian, then this result corresponds to the coupling 
term: 
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f{uj)S{Lu) ^^Jd^xJ d^x' fI{x)D^{x,x' ;w)F2{x')f{L 



(24) 



r 



Here Fi{x) = Q4o{x) and F2{x') = Q2oix'). For the 
pure hexadecupole giant resonance both Fi and F2 must 
be made equal to Q40. In general, mixed mode analysis 
produces strength functions with less pronounced peaks 
and can be used when computational time saving is nec- 
essary. 
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FIG. 7; The imaginary part of the response function corre- 
sponding to the isoscalar and isovector octupole, and isovector 
dipole are shown for ^^0 using the Sgll force. 



spherical RPA calculation using the same effective in- 
teractions The spherical calculation finds an ad- 
ditional broad peak for SkM* and Sgll centered at an 
energy of about 27 — 28 MeV. This peak is weaker than 
the 3 peaks at lower energies and is not observed in the 
three-dimensional linear response calculation. In Fig. 6 
one can see that the response is approximately purely 
negative, but that there are some violations of this fea- 
ture. 

The isovector octupole response is less clear than the 
other collective modes. It is also possible that resonances 
seen in the isoscalar octupole response may also appear in 
the isovector response due to couplings. In this case the 
strength of the peak is expected to be most prominent in 
its primary channel. This is most easily observed in Fig. 7 
for the Sgll force. In this case the lowest resonances at 
about 8 — 9 MeV and at 14 MeV are very close to the 
lowest peaks in the isoscalar case but with substantially 
reduced strength, whereas the most prominent isovector 
peak is at about 20 MeV. The sorting of the isoscalar 
versus isovector mode is more complicated for the SkM* 
force. 

The isovector dipole response resulting from the linear 
response calculations is not as prominent as the isoscalar 
responses. This is most likely due to the absence of a 
strong collective dipole resonance for this nucleus. For 
the SkM* case, there are some peaks centered around 
17 — 18 MeV, while for Sgll the peaks range about 18 — 
21 MeV. 



IV. SUMMARY AND CONCLUSIONS 



B. Octupole and Dipole Modes 

The octupole and dipole giant resonances are not sym- 
metric about the z = plane and hence these nodes do 
not couple with the symmetric quadrupole and hexade- 
cupole modes. The isovector octupole moment is defined 



130 



^(3) J 



(25) 



In Fig. 6 the isoscalar dipole, isovector octupole, and 
isoscalar octupole responses are shown for ^^O using the 
SkM* interaction. In Fig. 7 the same quantities are 
again plotted using the Sgll force. The isoscalar octupole 
mode shows three sharp resonance structures at about 
6-7 MeV, 13 - 14 MeV, and 18 MeV, for the SkM* 
force, whereas the Sgll resonances are at the slightly 
higher energies. These three peaks are also found in a 



A method for evaluating the linear response theory us- 
ing TDHF is formally developed and implemented. This 
method allows one to construct the dynamic calculation 
directly on top of the static Hartree-Fock calculation. 
Therefore, by performing a sophisticated and accurate 
three dimensional static Hartree-Fock calculation, we 
have a correspondingly accurate and consistent dynamic 
calculation. A coherent description of static ground state 
properties, such as binding energies and deformations is 
given along with a description of the collective modes of 
nuclei. 

A three-dimensional collocation Basis-Spline lattice 
representation is used, which allows for a much more ac- 
curate representation of the gradient operator and hence 
a correspondingly accurate overall calculation. Exam- 
ple calculations of two spherical systems (^®0, *°Ca) and 
of a system with prolate quadrupole deformation ('^^S) 
are presented for the response functions corresponding to 
various isoscalar and isovector multipole moments. The 
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SkM* case for the axial isoscalar quadnipolc mode gives 
excellent results, obeying expectations from the theory 
and also satisfying the energy weighted sum rule. 



bilt University. Some of the numerical calculations were 
carried out on IBM-SP supercomputers at the National 
Energy Research Scientific Computing Center (NERSC). 
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